function u0 = lagrange22(A,F,B,G,epsilon,tol,max_it)
%function u0 = lagrange22(A,F,B,G)
%u0 is the result of an iterated solution with matrices of reduced size
%u0=[A;B]\[F;G];
if nargin == 4
    epsilon = 1.0E-6;
    tol = 1.0E-10;
    max_it = 10;
end
% [p,n] = size(B);
cvg = 0;
it_count = 1;
% oldErr = inf;
% newErr = 10000; 
L  = 1/epsilon*(B'*B);
K  = A+L;
bg = 1/epsilon*(B'*G);
u0 = K\(F+bg);   % important
% condest(K)
% u0 = cgs(K,(F+bg),1e-4, 20);

% u1 = bicg(A+L,A*u0+bg,tol,500);

t0=cputime;
while ~cvg && (it_count <= max_it)
  u1 = K\(A*u0+bg);
%   u1 = cgs(K, (A*u0+bg), 1e-5, 500);
  %
%   oldErr = newErr;
  newErr = norm(u1-u0,inf)/norm(u1,inf);
%   disp(['The error between iterations is ', num2str(norm(u1-u0,inf)),
%   '.'])
  disp(['The error between inner linear iterations is ', num2str(newErr), '.'])
  cvg = newErr <=tol;
  it_count = it_count + 1;
  u0 = u1;
end;
t1 = cputime;
fprintf('Number of linear iterations is %d. Iteration Costs %3.20f.\n', it_count, t1- t0);
